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ABSTRACT 

We obtain analytical expressions for an effective potential of interaction 
between two- and three-dimensional (2D and 3D) solitons (including the case 
of 2D vortex solitons) belonging to two different modes which are incoherently 
coupled by cross-phase modulation. The derivation is based on calculation 
of the interaction term in the full Hamiltonian of the system. An essential 
peculiarity is that, in the 3D well as in the case of 2D solitons with 

unequal masses, the main contribution to the interaction potential originates 
from a vicinity of one or both solitons, similarly to what was recently found in 
the 2D and 3D single-mode systems, while in the case of identical 2D solitons, 
the dominating area covers all the space between the solitons. Unlike the 
single-mode systems, stable bound states of mutually orbiting solitons are 
possible in the bimodal system. 
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1 INTRODUCTION 

Interaction between self-focussing cylindrical beams (spatial solitons) in bulk 
nonlinear media is a problem of obvious interest both by itself and for applica- 
tions. This interaction was studied experimentally and by means of numerical 
simulations in photorefractive media [J, and was simulated in detail in an 
isotropic model with the second-harmonic-generating (SHG) nonlinearity 0. 
In the latter model, it was demonstrated that the spiraling bound state of 
two cylindrical beams is unstable. A general analytical expression for a po- 
tential of interaction between far separated two- and three-dimensional (2D 
and 3D) solitons was very recently derived in JJJ. The potential also predicts 
an instability of the orbiting bound state of two solitons. 

A very convenient model for the study of the multidimensional solitons 
and their interactions is the cubic- quintic nonlinear Schrodinger (CQNLS) 
equation, in which the cubic nonlinearity is self-focusing, giving rise to the 
beams (2D solitons) or "light bullets" || (3D solitons), while the quintic 
term is self-defocusing, precluding the wave collapse in the model: 

iu t + V 2 m + \u\ 2 u — g\u\ A u = 0. (1) 

The coefficients in ([I]), except for g, can be set equal to 1 by means of scale 
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transformations, whereas g is the quintic-to-cubic nonlinear susceptibilities 
ratio. In the application to nonlinear optical media, the temporal variable t 
in (|l|) must be replaced by the propagation distance z, while the role of the 
third transverse variable is played by the "reduced time", t — z/V gT , V gr being 
the mean group velocity of the carrier wave (this implies that the temporal 
dispersion in the medium is anomalous) || . Finally, the Hamiltonian of this 
model is 

H U = J {\Wu\ 2 - l -\u\^ 9 -\u\^dv. (2) 

Vortex beams, with the vorticity ("spin") s — 1, and interactions between 
them, described by Eq. ([!]), were simulated by Quiroga-Teixeiro and Michinel 
§ . A remarkable result is the numerically discovered robustness of the vortex 
beams (which were found to be strongly unstable against azimuthal pertur- 
bations in the SHG model QJ). Note that the model (|l|) is not merely the 
simplest one that gives rise to stable multidimensional solitons: according to 
experimental data |J, the combination of the focusing cubic and defocusing 
quintic terms adequately represents the nonlinear optical properties of some 
real materials. 

The effective potential of the intersoliton interaction derived in 0] applies 
to a wide class of models, including Eq. ([!]). However, it does not apply to 
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bimodal systems including two equations with incoherent nonlinear coupling 
between them (mediated by the cross-phase modulation in nonlinear optical 
media), in the case when the two solitons (beams) belong to different modes. 
The simplest bimodal generalization of the model based on the Hamiltonian 
(H) is furnished by the Hamiltonian H = H u + H v + H int , where H v is the 
same expression as (§) with the field u substituted by another field v, and 
the interaction part of the Hamiltonian is 
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J [-f3\u\ 2 \v\ 2 + a (\u\ 4 \v\ 2 + \u\ 2 \v\ 4 )] dr, (3) 



with, generally speaking, arbitrary positive constants a and (3. The full 
Hamiltonian of the bimodal system gives rise to the equations 

iut + V 2 u+ (\u\ 2 + [3 \v\ 2 ^j u - (g\u\ 4 + 2a\u\ 2 \v\ 2 + a\v\ 4 ^ju = 0, (4) 
iv t + V 2 v + (\v\ 2 + 13 |m| 2 ) u - [g\v\ A + 2a\u\ 2 \v\ 2 + a |u| 4 ) v = 0. (5) 

Commonly known examples of optical bimodal systems are provided by two 
orthogonal polarizations of light, or two light waves with different carrier 
wavelengths ||. In the latter well as in the former one with the 

circular polarizations, the cubic cross-phase-modulation coefficient is /3 — 2. 
In the case of two linear polarizations, (3 = 2/3 (and the usual assumption is 
to drop additional four- wave mixing terms). The constant a is left here to 



be arbitrary, but, in the most interesting cases, the second term in (Q) will 
only produce a small correction to the effective interaction potential. 

As well as in the case of the single-mode system, the interaction between 
solitons in different modes depends on the separation between them, but, 
unlike the single-mode case, it is not sensitive to a phase difference between 
the solitons, hence the interaction is expected to be simpler than inside the 
same mode. The objective of this work is to find an effective potential of 
the interaction between 2D and 3D solitons in the bimodal system, including 
the cases when the interacting solitons are both identical and different (the 
interaction between 2D solitons with different vorticities is also included). 
The interaction between identical 2D beams was recently considered in [TO 



but using an artificially introduced Gaussian ansatz for the soliton. As well 
as it was done in [[| for the single-mode system, in this work we find the 
interaction potential in a general analytical form. However, the derivation is 
essentially different from that developed in 0]; in particular, the derivation 
proves to be very different for the 2D and 3D cases, while in the single-mode 
system these two cases were very similar. In section 2, we derive the potential 
for the interaction between 2D solitons (spatial beams) with unequal masses. 
We explicitly consider two limit cases, when the masses of the interacting 



6 



solitons are very different or nearly equal. The latter case demonstrates 
a singularity in the limit of equal masses, therefore the interaction between 
identical 2D solitons should be considered separately, which is done in section 
3. The result, and the way to obtain it, turn out to be drastically different 
from the case of unequal masses: when the masses are not equal, a dominating 
contribution to the effective interaction potential is produced by a vicinity 
of the soliton having a larger mass (which is similar to the situation for the 
single-mode system ||), while, in the equal- mass case, a dominating area 
is located between the solitons, in contrast with the case of the single-mode 
system. In section 4, the potential is derived for the 3D solitons with equal 
masses. In this case, the derivation is similar to that for the 2D solitons with 
the unequal masses, but it gives rise to an additional logarithmic factor. 

The results are summarized in section 5, where we conclude, in particu- 
lar, that the obtained potentials give rise to two bound states of the solitons 
orbiting around each other, one of which is stable (which was already con- 
cluded in []10[), on the contrary to the orbiting states in the single-mode 
models Q (including the SHG one 0), which are all unstable. This differ- 
ence, which is, obviously, very important for applications, is due to the fact 
that, in the bimodal system, the interaction potential does not depend on 
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the phase difference between the two solitons. 



2 THE INTERACTION BETWEEN DIFFER- 
ENT TWO-DIMENSIONAL SOLITONS 

A general 2D stationary solution to Eq. (Q), with an internal frequency 
<jj = —fi 2 , is looked for in the form 

u s = exp (ifi 2 t + isO^j U (r), (6) 

where the s is an integer spin (vorticity), and a real function U{r) satisfies 
the equation U" + r~ 1 U' — s 2 r~ 2 + U 3 — gU 5 = fi 2 U, which can be easily 
solved numerically ||. A soliton solution is defined by its asymptotic form 
at r — ► oo, 

U(r) w A s ([i) r" 1/2 exp(-/ir), (7) 

where the amplitude A s (fi) is to be found numerically. We will consider a 
situation with the solitons of the form ([]), (|7|) in each mode u and v, that 
have, generally, different spins s\ and S2 and different frequency parameters 
fix and fi2 which determine effective masses of the solitons. A size of the 
soliton can be estimated, pursuant to Eq. (0), as fi^ 1 . We will consider 
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the case when the separation between the solitons is much larger than their 
proper sizes, i.e., Rfii t 2 ^> 1- 

The interaction Hamiltonian (D) allows one to define an effective inter- 
action potential for two separated solitons, approximating the two-soliton 
configuration by a linear superposition of two isolated solitons, and sub- 
stituting it into (|D [p| . This approach requires actual calculation of the 
integrals in (||), which can be done for 2D solitons in an exact form only 
in exceptional cases (see, e.g., [0), another drawback being that a distor- 
tion of the "tail" of each soliton due to its interaction with the "body" of 
the other one is ignored. In the work [TO] , the necessary integral was eval- 
uated, assuming a Gaussian ansatz for the isolated solitons. However, the 
corresponding effective interaction potential, decaying ~ exp (—const • R 2 ), 
was actually produced by the ansatz rather than by the model. In fact, the 
potential must decay as exp (—const • R), see below. 

In the work |4[], another approach to the calculation of the effective poten- 
tial was developed for the single-mode systems, following the way elaborated 



earlier for ID solitons in [13]. This method did not require knowing inter- 
nal structure of the soliton, and did not imply neglecting the distortion of 
each soliton's tail due to its overlapping with the other soliton. All that is 
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necessary to know about the individual solitons for the application of this 
method, is only the asymptotic amplitudes A s (fj) in (|7|). Here, we will apply 
a similar method to the bimodal system @, (H), although technical details 
will be essentially different from those in the case of the single-mode systems. 

We start, still assuming the linear superposition of the two solitons u Sl 
and v S2 (see Eqs. ([]) and (0)) with widely separated centers, and setting, for 
the definiteness, /ii > \i 2 . Because of the exponential decay of the fields, it 
is obvious that a dominant contribution to the interaction potential (0) will 
be produced by a vicinity of the soliton with fi = 

First, we consider the case /12 <C /ii, i.e., a light soliton (beam) interact- 
ing with a heavy one. Substituting the field v for the light soliton by the 
asymptotic expression ([7p, and neglecting its small variation over the size of 
the narrow heavy soliton, we can easily cast the expression for the interaction 
potential into the form 

H int « ^ 2 ( /i2 ) J R- 1 exp(-2 /U2 i?) 

= (-Pm 1 + am 1 )A 2 S2 (fi 2 )R~ 1 exp(-2/2 2 R) , (8) 

where mi = J \u Sl (r;fii)\ 2 dr and mi = / |u Sl (r;/ii)| 4 dr are two integral char- 
acteristics of the heavy soliton, mi being, in fact, its effective mass. Thus, 
both attraction and repulsion between the light and heavy solitons may take 
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J \u Sl (r;ni)\ 2 dr + a J |u Sl (r;/ii)| 4 dr 



place, depending on the sign in front of the expression 



Another interesting case is fi\ — /12 = A/i <C ^2 = A* (i- e -; t ne interac- 
tion between nearly identical solitons, provided that s% = 82). In this case, 
following |4[], we assume that, in terms of the polar coordinates (r, 9) with 
the origin at the center of the heavier (first) soliton, a main contribution to 
H int comes from the distances <C r <C R, where both solitons may be 
approximated by the asymptotic expressions (0). Then, it is straightforward 
to obtain the following expression corresponding to the first term in (Q) (cf. 
the corresponding expressions and Fig. 1 in [|J): 



Here, the small difference A/z, and the difference between R and the exact 
distance from an integration point (r, 9) to the center of the second (lighter) 



soliton, \J (R + r cos#) 2 + r 2 sin 2 9 , are neglected everywhere, except for the 

argument of the exponential function. Making use of the expansion 








(10) 



and of the formula J 2?r exp (— 2/ir cos 9) d9 ~ 



Tr/firexp (2 fir), valid for fi r ^> 1, 
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we can simplify the integral @ to the form 

H int « -/5<( /U )^ 2 ( /U ) v /^Vi?- 1 exp(-2 /U i?) /' 00 r- 1 / 2 exp[-2(A/i)r] C /r 

w 

= -7T/5(2 /U A M )- 1/2 <( /i )^ 2 ( / i)i?~ 1 exp(-2 /Uj R). (11) 

For small A/x, the integral in (|TTD is dominated by a contribution from the 
region r ~ 1/A/i ^> which justifies the use of the asymptotic approxi- 
mation (0) for the field u(r). A contribution from the second term in the full 
expression (^), evaluated in the same approximation, demonstrates the same 
dependence on the separation R, but without the large multiplier (A/i) -1 / 2 , 
therefore this only a small correction to (|TT1) . 

3 ATTRACTION BETWEEN IDENTICAL 
TWO-DIMENSIONAL SOLITONS 

The expression (O) diverges in the most interesting case A/x = 0, which 
corresponds to the interaction between identical solitons (provided that si = 
S2 = s). The divergence suggests that a region dominating the interaction 
potential is not that around the soliton, as was the case both in the previous 
section for the case of A/i ^ 0, and, for the identical solitons, in the single- 
mode system Q, but a wider region between the two solitons. To calculate 
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U int in this case, we use the Cartesian coordinates (x, y) defined so that the 
centers of the two solitons are placed at the points (±i?/2,0). Then, using 
once again the asymptotic expressions (0) for both u- and f-solitons with 
fi 1 = fi 2 = /i, the interaction potential corresponding to the first term in Eq. 
(0 ) is given, after obvious transformations, by 

-1/2 



exp 



(e 2 + r/ 2 + l/4) -e 



-2^R (V (£ + 1/2Y + r ? + y/(Z- 1/2) 2 + rf 



(12) 



where £ = x/R, and r\ = y/R. Because the parameter /ii? is large ac- 
cording to the underlying assumption, the integral ([12]) is dominated by a 



contribution from a vicinity of points where the argument of the exponential 
function has a maximum. An elementary analysis shows that the maximum 
is attained not at isolated points, but rather at the whole segment |£| < 1/2, 
i] = 0. Expanding the integral in small t] 2 in a vicinity of this segment, we 
approximate Eq. ([H|) by an integral that can be easily calculated: 



r + l/2 /I \-l r+oo ( 1 

#int = -(3A 4 S J ^ 2 d£ -e) J ^ d V exp -2fiR \l + - 



Tf 



V 21/4-f 



(13) 



Comparing this result with that (|il|) for the solitons with different masses, 



we conclude that the divergence in the latter expression at A/x — > im- 
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plies replacement of the pre-exponential factor R^ 1 by a larger one, R~ x l 2 . 
Evaluating the second term in (|3|) in the same approximation, we conclude 
that it yields an expression differing from (|i~3|) just by the factor R~ l instead 
of i? -1 / 2 , i.e., a small correction to ([131) . Note that, unlike the interaction 
between heavy and light solitons, which may have either sign, the nearly 
identical or identical solitons always attract each other, cf. Eqs. (||), (|i~TD, 
and flTJ). 

4 ATTRACTION BETWEEN THREE-DIMENSIONAL 
SOLITONS 

In the 3D case, we consider only the solitons without the internal "spin", 
i.e., with s = 0. The 3D soliton has the form u = exp(z/i 2 t) a(r), with the 
asymptotic form a(r) « Ar~ l exp(— fir) at r — > oo, cf. Eqs. (||) and ©• 
Substitution of this asymptotic expression for the fields u and v into the 
integral in the first term of Eq. (||) around each soliton (cf. Eq. (|| )) and 
using the expansion ([[(]) yield 

/•OO /*7T 

H int = -Anf3A A R- 2 / r 2 dr / sin d6 r~ 2 exp [-2/xr -2/i(R + r cos 9)] , 
Jo Jo 



14 



where an extra factor 2 takes into regard the fact that, in the case of identical 
solitons, one has equal contributions from the vicinity of both solitons. After 
elementary integration over 9, we arrive at an expression containing a formal 
logarithmic singularity, 

POO 

flint = -27r/?AV 1 ^" 2 exp(-2/i J R) / r- l dr. (14) 

Jo 

Actually, the lower and upper limits of the integration are, respectively, 
~ and ~ R, so that, with the logarithmic accuracy, the final expression 
for the effective interaction potential in the 3D case becomes (cf. ([3~3|)) 

fl int = -27r/5A 4 /i" 1 J R" 2 exp (-2/xfl) ln(/xfl). (15) 

As for the contribution from the second term in (131), it has the same depen- 



dence on R as (]T5|) , but without the large logarithmic factor, so that in this 
case too, it is a small correction only. 

The above analysis was done for identical solitons. If the solitons have a 
small mass difference, corresponding to a small difference A/z, the interaction 
potential is given by essentially the same expression (|15"D, except for a factor 
of 2, which is absent if A/i • R > 1, when the calculation of H int is dominated 
by a vicinity of one soliton only. 
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5 CONCLUSION 

In this work, we have derived analytical expressions for an effective potential 
of interaction between two- and three-dimensional solitons (including the 
case of the two-dimensional vortex solitons) belonging to two different modes 
which are incoherently coupled through cross-phase modulation in models of 
media with the self-focusing cubic and self-defocusing quintic nonlinearities. 
The derivation was based on the calculation of the interaction term in the 
full Hamiltonian of the system. An essential peculiarity is that, in the 3D 
case, as well as in the case of 2D solitons with unequal masses, the main 
contribution to the interaction potential originates from a vicinity of one or 
both solitons, similarly to what was recently found in the 2D and 3D single- 
mode systems Q , while in the case of identical 2D solitons, the dominating 
area covers all the space between the solitons. Except for the case of the 
interaction between light and heavy solitons, which may have either sign, 
the solitons always attract each other. 

The attraction between the solitons may give rise to their orbiting bound 
states in the 2D and 3D cases (in the latter case, it is assumed that the two 
solitons move in one plane). Orbiting of incoherently interacting 2D solitons 
was experimentally observed in a photorefractive medium [TJ. Numerical 
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simulations and analytical arguments presented in || and 0] demonstrate 
that the orbiting bound states of the 2D solitons in the single-mode systems, 
including the SHG model, are unstable. However, it was recently pointed 



out in |T0[ that they might be stable in the bimodal system. Indeed, for 
the orbiting state, the interaction potential (|TT|) , flT3|) , or (|15D must be sup- 
plemented by the centrifugal energy E ci = (M 2 /2m) R- 2 , where M is the 
angular momentum of the soliton pair, and m is the soliton's effective mass. 
A straightforward consideration of the net potential, H irit + E c f, demonstrates 
that it may have two stationary points, the one corresponding to a smaller 
value of R being a potential minimum that gives rise to a stable orbiting 
state. The instability of similar states in the single-mode systems is due to 
the fact that, in those systems, the interaction potential also depends on the 
phase difference between the solitons, the effective mass corresponding to 



the phase-difference degree of freedom being negative This, eventually, 
made the existence of stable stationary points of the effective interaction 
potential impossible. 

This work has been supported by INTAS under grant No: 96-0339 
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